THE CELLULAR DYNAMICS OF BONE REMODELING: 
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Abstract. The mechanical properties of vertebrate bone are largely determined by a process which involves the 
complex interplay of three different cell types. This process is called bone remodeling, and occurs asynchronously 
at multiple sites in the mature skeleton. The cells involved are bone resorbing osteoclasts, bone matrix producing 
osteoblasts and mechanosensing osteocytes. These cells communicate with each other by means of autocrine and 
paracrine signaling factors and operate in complex entities, the so-called bone multicellular units (BMU). To investi- 
gate the BMU dynamics in silico, we develop a novel mathematical model resulting in a system of nonlinear partial 
differential equations with time delays. The model describes the osteoblast and osteoclast populations together with 
the dynamics of the key messenger molecule RANKL and its decoy receptor OPG. Scaling theory is used to address 
parameter sensitivity and predict the emergence of pathological remodeling regimes. The model is studied numeri- 
cally in one and two space dimensions using finite difference schemes in space and explicit delay equation solvers in 
time. The computational results are in agreement with in vivo observations and provide new insights into the role 
of the RANKL/OPG pathway in the spatial regulation of bone remodeling. 
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1. Introduction. The vertebrate skeleton plays a crucial role in providing mechanical sup- 
port as well as a ready source of calcium and other important minerals. Physical loading of the 
skeleton causes stresses which can lead to local micro-damage in the bone tissue. Similarly, if the 
calcium level in the blood drops below a certain threshold, systemic regulators such as hormones 
transmit the order to release calcium through removal (resorption) of bone tissue. In both cases, the 
resorbed spaces have to be filled with sound tissue in order to restore the structural integrity. This 
joint process of bone destruction and re-growth is referred to as bone remodeling, and is realized by 
complex multicellular entities, the so-called bone multicellular units (BMU). Each BMU consists of 
several interacting cell types and a whole variety of biochemical signaling factors. The importance 
of remodeling becomes apparant when considering the implications of its malfunctioning. Deficient 
or even absent remodeling of micro-damage can lead to macroscopic bone fractures, and pathologies 
in the BMU functioning are largely responsible for diseases such as osteoporosis and rheumatoid 
arthritis 25 . 

The various physiological and pathological aspects of BMUs have been studied by both exper- 
imentalists and clinicians for well over 40 years f29^. However, due to a general lack of conclusive 
in vivo experiments - so far mainly consisting of histological sections of dead bone tissue - several 
phenomena remain poorly understood. The difficulty and costs for in vivo experiments suggest 
that there is great potential for mathematical modeling in this field. So far, several research groups 
have modeled the local strain fields in bones [331 131] as well as the temporal sequence of local bone 
destruction and re-growth at the cellular level [23l[24l[22. In essence, the latter models success- 
fully capture the local bone cell dynamics in physiological settings and are even able to describe 
certain pathologies. However, the functioning of a remodeling unit strongly depends on its spatial 
organization and therefore, a purely temporal model cannot provide a complete description of the 
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BMU. To address this, we develop here a novel spatio-temporal model of a single remodeling unit, 
describing the dynamics of both the involved bone cell populations as well as the relevant signaling 
pathways. The model consists of five nonlinear partial differential equations (PDE) and is based 
on a continuum assumption for the cell populations. 

In section [2] we first give an outline of the relevant biology, thereby focusing on the three types 
of bone cells (osteoclasts, osteoblasts, osteocytes) and the most important biochemical factors (the 
RANK/RANKL/OPG pathway). Once these concepts are established, we begin the model devel- 
opment in section |3] by introducing a previous temporal model by Komarova et al. [23l [24] . Given 
the complexity of the underlying biological system - involving endocrine signaling, cell motion, 
fluid diffusion etc. - some simplifying assumptions are necessary in order to develop a compact 
and closed spatio-temporal model. The model is developed in an abstract setting independent of 
the spatial dimension, but can be applied to one (ID), two (2D) or three (3D) dimensions. In 
section |4] we present the ID case, use scaling theory to gain insight into parameter sensitivity, and 
present experiments focusing on the different pathological regimes. The biologically more relevant 
2D case is then discussed in section [5] and a selection of two physiological remodeling experiments 
is presented. The results of the 2D experiments provide a model validation as well as new insights 
into the role of the RANK/RANKL/OPG pathway in the spatial regulation of bone remodeling. 

2. The biology of bone remodeling. Bone remodeling refers to the combination of bone 
destruction and subsequent re-growth. It is a coordinated process of three different cell types that 
interact by means of several biochemical factors. Furthermore, mechanical strains play an important 
role in the stimulation and steering of remodeling units. The following outline is focused on the 
model-relevant mechanisms and we refer to [29l [3T] for detailed reviews. 

2.1. The bone cells. Three different cell types are involved in remodeling. 

• Osteoclasts jJO, i4j are cells which resorb mineralized bone tissue while moving along the bone 
surface. They are formed by cell differentiation from stem cells in the bone marrow and have a 
life span of roughly 10 days. A key stimulator for osteoclast differentiation and activation is a 
molecule called RANKL (the receptor activator of nuclear factor kB). 

• Osteoblasts [15] are cells which fill the previously resorbed trench with osteoid, the organic 
part of the bone tissue. Later on, osteoid mineralizes and the remodeling process is complete. 
Osteoblasts differentiate from stem cells in the bone marrow, they do not move along the bone 
surface, and they express the messenger molecule RANKL and its decoy receptor OPG (osteo- 
protegerin). After approximately two weeks, osteoblasts either die or differentiate into osteocytes 
and get buried alive in the new bone tissue. 

• Osteocytes [T31 [3] differentiate from active osteoblasts and are connected with each other to 
form a large network of active cells within the bone tissue. This network is believed to propagate 
information, to localize damage sites and micro-strains, and to play an important role in the 
process of mechanotransduction. 

The three cell types communicate by means of autocrine signaling (communication among cells of 
the same type) and paracrine signaling (communication among cells of different types) . Generally, 
the bone cells and their messengers operate locally in well-confined remodeling units, the BMUs. 
These units operate for up to 12 months in a row, thereby by far exceeding the individual cell's life 
spans. The progression of a BMU across the bone can be summarized as follows: 
Step 1) Initially, 10—20 osteoclasts are recruited to the initiation site and resorb the old bone tissue. 

Once the tissue is removed, the osteoclasts move on and keep on resorbing while traveling 
at a speed of 20 — 40/im per day [iniHn]- During the whole remodeling process, they stay 
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together in a spatially well-confined aggregation (cutting cone). Dead cells are continually 
replaced by new ones so that the population size remains approximately constant. 
Step 2) Once the osteoclasts have resorbed the bone tissue, they recruit 1000-2000 osteoblasts that 
fill the previously resorbed trench with new bone matrix (closing cone). Osteoblasts are 
much less efficient than osteoclasts and the bone formation takes roughly 10 times longer 
than the resorption. 

Step 3) Finally, the new bone matrix mineralizes and osteoblasts either die or differentiate into 
osteocytes. 

There are two kinds of bone tissues. Cortical tissue is dense and compact and forms the outer 
surface of bones. Trabecular tissue fills the inner cavity with a honeycomb-like structure, consisting 
of irregularly shaped spicules (trabeculae) endowed in bone marrow. Remodeling takes place in both 
cortical and trabecular bone and the difference in the respective BMU progressions is geometrical 
rather than biological in nature: whereas the BMU has to dig a complete tunnel to penetrate the 
compact cortical tissue, it can move along the surface of the trabeculae, thereby only digging a 
half-trench. Figure [23] illustrates the temporal sequence of the remodeling steps on a trabecula. 




Fig. 2.1. A schematic, not-to-scale representation of a BMU moving along a micro-fracture on a piece of 
trabecular bone. Osteoclasts resorb the bone in form of a cutting cone and osteoblasts subsequently fill the resorbed 
space with new bone matrix. Bone cells interact by means of cytokines and growth factors and osteoblasts differentiate 
into osteocytes. 



2.2. The biochemical factors involved in remodeling. The coordination of osteoclasts, 
osteoblasts and osteocytes within a BMU is realized through a sophisticated communication system, 
which consists of various autocrine and paracrine signaling pathways involving numerous coupled 
effectors. However, the multiple actions attributed to some of these effectors make it hard to 
identify the actual key players and to predict the cumulative dynamics of the coupling. Figure 



2.2 summarizes the major control pathways in the remodeling process and identifies the respective 



messenger molecules. Among the multiple messengers involved, RANKL and OPG have been shown 
to play critical roles in both physiological bone remodeling and in the development of diseases 
PSI ini m] . RANKL is a cytokine produced in either membrane-bound or soluble form by cells 
of the osteoblast lineage, prominently by osteocytes and osteoblasts. Several studies have shown 
that RANKL is up-regulated in situations associated with increased bone remodeling, such as PTH 
treatment [17], mechanical stimulation [20], as well as fractures [18]. RANKL binds to RANK 
receptors on the surface of osteoclastic cells and has a stimulatory impact on the differentiation 
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of osteoclast precursors and the subsequent activation of mature osteoclasts into active, resorbing 
cells. On the other hand, the molecule OPG, produced by mature osteoblasts [S], acts as a 
decoy receptor of RANKL, i.e. it inhibits RANKL by forming RANKL-OPG complexes. Since the 
presence of OPG means less RANKL-RANK binding and hence less osteoclast stimulation, a high 
RANKL/OPG ratio favors bone resorption whereas a low ratio down-regulates osteoclastic activity. 
The RANK/RANKL/OPG pathway is also known to be employed by systemic regulators such as 
parathyroid hormone (PTH) and vitamin D to regulate the resorption activity. Note finally that 
the spatial separation of the different RANKL and OPG sources indicates that in addition to the 
local ratio of the chemicals, their spatial distribution plays an important role, too. 



Fig. 2.2. Cells and biochemical factors known to play a role in the remodeling process of bone. The cells are 
osteoclasts (OC), osteoblasts (OB), osteocytes (OCY), and their respective precursor cells. Solid lines stand for 
positively balanced processes (cell differentiation and production of chemicals/tissue) and dotted lines for positively 
balanced regulations (autocrine/paracrine stimulation). The (-) next to an arrow indicates a negatively balanced 
process or regulation. 

2.3. The mechanical effects: microscopic strains and fractures. There are two differ- 
ent remodeling modes, targeted and random remodeling. Whereas the former mode aims at damage 
removal by means of local micro-fracture reparation, the latter serves the purpose of damage pre- 
vention: old - but not necessarily damaged - tissue is continually renewed across the skeleton 
to prevent fatigue damage. Both remodeling types rely on steering mechanisms that ensure that 
BMUs are guided towards damage sites and move in a way that minimizes structural instability due 
to ongoing bone erosion. The concept of targeted steering is based on established evidence that the 
presence of micro-fractures leads to creation of new BMUs and attraction of already existing BMUs 
[SlITj- On the other hand, it has been suggested that strain-derived canalicular fluid flow is responsi- 
ble for osteoclast activity and motility in the cutting cone of the BMU [B], leading to strain-derived 
steering. In particular, this steering mechanism ensures that BMUs move along the principal strain 
axis of the bone and hence optimize its robustness at anytime of the remodeling process. Both steer- 
ing mechanisms rely on mechanical features that need to be translated into cell signals to attract 
BMUs. Recent investigations show that there is a unifying mechanism of mechanotransduction for 
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both damage and strain, mediated by osteocytes. In fact, both mcchanicany damaged osteocytes 
and osteocytes exposed to fluid shear stress have been shown to express RANKL [2 [37] . Since 
RANKL is a potent osteoclast stimulator, this allows mechanically stimulated osteocytes to attract 
BMUs and hence guide them towards damage sites and along the principal stress directions. 

3. The Mathematical Model. In this section, we develop a mathematical model describing 
the spatio-temporal evolution of a single BMU at cellular level. The overall goals of this model are 
the following: 

• To describe the distinctive spatial and temporal features of the cutting cone and the BMU 
movement. 

• To link the key biochemical factors RANKL and OPG with the known population dynamics 
of bone cells. 

• To test the model on experimental findings and suggest new experimental studies. 

Since we develop a model that can be considered in one, two and three space dimensions, we do 
not specify its dimension explicitly and denote it simply by n, where n = 1,2,3. The ID and 2D 
versions of the model presented in this article are particularly suited for the description of trabecular 
remodeling, and the restrictions of their applicability to cortical bone will be discussed in section |6] 
The major modeling assumptions can be summarized as follows: 

• We focus on trabecular remodeling, more precisely on the dynamics of a BMU moving 
across a single trabecula. 

• The trabecula is locally flat enough so that we can neglect curvature. 

• We make a continuum assumption for the cell population densities, i.e., we shall be modeling 
cell densities rather than individual cells. 

• The BMU evolves along the surface of the trabecula and the depth of the resorbed trench 

lOfim) is small in comparison to its width {'^ 500/xm). 

• Of the several cell types involved in remodeling - osteoblasts, osteoclasts, osteocytes and 
their respective precursors - we only consider osteoblasts and osteoclasts as state- variables. 

• The trabecula is endowed in bone marrow which can be considered as a reservoir of pre- 
cursor cells. 

• Among the multitude of biochemical factors, only RANKL and OPG are modeled explicitly, 
the rest of the factors - such as TGF-/3, IGF, M-CSF and nitric oxide - are captured in 
nonlinear interactions. 

• The canopy of bone lining cells separating the BMU from the bone marrow [12] ensures 
that the loss of chemicals by vertical diff'usion is negligible. 

• We model the elimination of OPG and RANKL through their mutual interaction only and 
do not include their natural decay rates. 

• The mechanical factors responsible for the BMU steering - microscopic strains and damages 
- are modeled implicitly in form of an appropriate RANKL distribution in the initial field. 
For the sake of simplicity, we will from now on refer to these distributions as micro- fractures, 
even though they might be caused by local micro-strains, see section [2.31 

Due to the complexity of the model, we proceed in three steps, starting off with a brief review of 
the temporal model introduced by |23|, I24j. In a second step, we introduce the spatial extension 
of the model as well as the RANKL and OPG fields. In a third step, we complete the model by 
adding appropriate initial and boundary conditions. 

3.1. Prior work: temporal model. The model suggested by Komarova et al. [231 El] is a 

temporal model describing the population dynamics of bone cells at a single point within the BMU. 
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Denoting the number of osteoclasts and osteoblasts by ui and U2 respectively, the cell dynamics 
are modeled as 



dtUi — aiu1^^U2^^ — f3iui 

dtU2 = a2uf'^uf^ - I32U2 



where and Pi are activities of cell production and death and all have units [day~^]. The four 
dimensionless parameters gij represent the effectiveness of the autocrine and paracrine interactions 
between the constituent cells. Let us now briefly discuss the various signaling factors gij, thereby 
making some restrictions appropriate to the spatio-temporal model we are finally aiming for. The 
factor gii represents the effectiveness of the autocrine interactions between osteoclasts and has been 
shown to control the overall remodeling dynamics [23]. Osteoclast-derived paracrine regulation 
of osteoblasts ((712) is the crucial link in the BMU coupling and its inhibition leads to negatively 
balanced remodeling |23j . Regarding the autocrine stimulation of osteoblasts (522), it is known that 
the latter express auto-stimulatory factors such as insulin-like growth factors IGF [9|. However, 
these factors do not influence the dynamical behavior of the BMU [24] and we assume here that 
they are negligible in comparison to the impact of gi2, i.e we set 522 = 0- Finally, osteoblast-derived 
paracrine regulation of osteoclasts is dominated by the RANK/RANKL/OPG pathway [5S1[55] and 
therefore the factor 321 plays an important role in the temporal model. However, since we will 
eventually develop a model that includes the RANKL and OPG flelds explicitly as state-variables. 



we can set 1721 ~ 0. After these simpliflcations, the system (3.1 1 reduces to 

dtU2 = a2uf'^ - /32U2- ' 



For gii < 1, the unique non-trivial flxed point {ui^ss,U2^ss) > (OjO) of equation (3.2) is a stable 
node. It is assumed that cells below the steady-state values Ui^ss are precursor cells which are less 
differentiated. In other words, they are not actively involved in the resorption and production of 
bone matrix, but participate in autocrine and paracrine signaling. Increases in Ui above Ui^^s are 
regarded as proliferation and differentiation of precursors into mature osteoclasts and osteoblasts 
that participate actively in the remodeling process. In this sense, the initiation of remodeling can 
be induced manually by increasing the number of osteoclasts above the equilibrium value, i.e. by 
choosing initial conditions ui(io) > ui^ss- Note that 1*2(^0) = U2,ss is sufficient because it is assured 
that osteoblasts are recruited by active osteoclasts. For all the subsequent numerical experiments 
we will choose the parameter gu < 1 such that {ui^ss,U2,ss) corresponds to a stable steady-state 
solution of ( |3.2| ). Together with the initiation procedure explained above, this implies that (wi , U2) > 
(ui^ss,U2,ss) for all t > to and hence we can ensure that the populations of active cells, denoted 
hereinafter by yi = Ui — Ui^ss, remain non-negative. Using the decomposition u,; = Ui^ss + yi, we 



can see that the system (3.2) actually describes the evolution of the active cell populations coupled 



to the constant precursor populations 

dtyi = ai{ui^ss + 1/1)^'-'- - Piiui^ss + yi), 



(3 3) 

dty2 = a2(wi,ss + yi)^'" - /32(W2,SS + y2)- ^ 



Even though our main interest is the evolution of the active cells in (3.3 1, we will henceforth use 



the equivalent version (3.2) for its more compact notation. The active cell populations are then 
easily recovered by subtracting the corresponding precursor populations Ui^^s from the solutions Ui 
of (|0). 
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3.2. The spatial extension. We use now the temporal model (3.1 ) as the basis for the spatial 
extension. The model developments in this section are independent of the spatial dimension and 
we avoid a specific choice by denoting all differential operators by their multidimensional symbols 
such as V and A. Later on we will discuss the ID case in section |4] and the 2D case in section [5] 
The units of the parameters introduced below can all be found in the Appendices |B] and [C] 

First, we switch to space-dependent state variables Ui{t) i— )■ Ui{x,t), where x G H C M" and 
the domain 11 is chosen large enough to avoid interactions with the boundaries (n — 1,2,3 is the 
spatial dimension). Note in particular that the Ui have now the units of a density [mm~"]. At 
the same time we introduce the RANKL and OPG fields as new state variables. They are denoted 
by 4>R{x,t) and (j)o{x,t) and have the units of a concentration [mol mm""]. To build up the 
final model we proceed now in two steps. First, we assume that the RANKL and OPG fields are 
known and analyze their impact on osteoclasts and osteoblasts. In a second step we introduce then 
the equations governing the spatio-temporal evolution of the RANKL and OPG fields themselves. 
Finally, we would like to emphasize that throughout the spatial extension the quantities Ui^ss refer 



to the steady-state densities of the temporal equation (3.2) and not to the steady-state solutions 
of the spatial equations. 

3.2.1. The impact of RANKL and OPG on osteoclasts and osteoblasts. RANKL is 
known to have an important impact on osteoclasts: it promotes their differentiation and activation 
and contributes together with other signaling molecules to the navigation (chemotaxis) of active 
cells [H [T^. On the other hand, the only impact of OPG on osteoclasts is indirect by means 



of RANKL inhibition. Accordingly, the osteoclast equation in (3.11 has to be augmented by two 
contributions only: 

dtui = aiuf " - /3iui - C V • {yiV^rd + h 9iyi) u, . (3.4) 
. ' ^ + <Pr 



ci 



C2 



The term CI describes the motion of active osteoclasts along the gradient of the RANKL field 
and C indicates the effectiveness of migration. The second term C2 represents the stimulating 
action of RANKL on osteoclasts via RANK-RANKL binding (fci is the corresponding reaction 
rate). This comprises both the differentiation of precursor cells into active osteoclasts as well as 
the steadily occurring renewal of nuclei in already resorbing cells [29 . We assume that the RANK 
receptors have a saturation threshold, hence the sigmoid function with A as the concentration of 
half-saturation. The Heaviside function 0{yi), defined as {d{x) = if x < 0, d{x) = 1 if x > }, 
ensures that stimulation takes place only in presence of active osteoclasts (yi), i.e. only osteoclasts 
(ui) in the cutting cone area are stimulated by RANKL. It is easy to verify that if ui(to) > "Ui.ss, 
then ui > ui^ss for all t > to, i.e. the population of active osteoclasts stays non-negative. Therefore, 



the same comments as in section 3.1 apply and equation (3.4) can, similarly to equations (3.2) and 
(3.3), be rewritten as an evolution equation for yi. 

Regarding osteoblasts, we assume that they are recruited by osteoclasts and do not move by 
themselves. Since RANKL and OPG have no significant impact on their dynamics, the U2 equation 



in (3.1) remains unaltered. 



3.2.2. Dynamics of RANKL and OPG fields. The evolution of the RANKL concentration 
4>ii is governed by production, diffusion and reaction. More precisely, RANKL is expressed by active 
osteoblasts, it spreads across the trabecula through diffusion and it binds to OPG as well as RANK 
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receptors on osteoclasts. In mathematical terms, the rate of change in time reads 



dt(t)R = an y2^tn+ krM4''r 



9{yi) ui - k'icjiRcjyo ■ 



(3.5) 



C3 



C4 



C6 



C5 



The RANKL source by active osteoblasts C3 is justified as follows: after the differentiation of 
precursors into mature osteoblasts, it takes a certain time until the cells start to produce RANKL 
[14j[34l[2]- The number of active osteoblasts at time t that are of age > tor older is e~^^*'^y2{x, t— 
tfi) and after absorbing the constant prefactor into the proportionality constant qr we obtain C3, 
where 2/2, tH = y2{x,t — tn). The second contribution C4 takes care of the porous diffusion which 
can vary between very low for membrane-bound RANKL and high for soluble RANKL. kj^ is the 
diffusion constant and the dimensionless exponent e_R > 1 reflects the porosity of the medium 
surrounding the BMU. Note that ii cr > 1, then an initially compactly supported RANKL field 
will stay compactly supported over time; this is not the case for the regular diffusion equation which 
is known to have infinite propagation speed. Since the BMU environment is very irregular and since 
the spreading cytokines are in steady interaction with the various constituents of the bone matrix as 
well as adjacent bone cells, the porous version with > 1 seems to provide a more plausible model 
for the RANKL field than the regular version with = 1. For a more detailed discussion of porous 
medium equations we refer to |12( lllj. The contribution C5 is due to the receptor-ligand binding 



of RANK and RANKL and is almost identical to C2 in equation (3.4), except for the different rate 
constant fc2- Note that ^2 contains information about several factors such as receptor density on 
osteoclasts and reversibility of the RANK-RANKL binding. Finally, the reaction term C6 models 
the RANKL-OPG binding with rate constant k^. 

Similarly to </)/{, the rate of change in the OPG field (po is also governed by the contributions 
of source, diffusion and reaction: 




(3.6) 



Similarly to C3 in equation (3.51, OPG is produced by mature osteoblasts with a time delay to 
such that to > tR [Ill[34, 2 . The contribution C8 for porous diffusion (eo > 1) is analog to C4 and 
the OPG-RANKL binding C9 is identical to C6. Note that the diffusion parameters of RANKL 
{f^Ry ^r) a-nd OPG {ko, eo) a-re not necessarily equal. In a physiological setting, RANKL is mainly 
membrane bound whereas OPG is soluble. 

3.3. The complete modeL Together with the evolution of the bone density z{x, t) - dimin- 
ished by active osteoclasts and augmented by active osteoblasts - equations (3.1), (3.4|, (|3.5[) and 



(3.6) yield the following nonlinear, time-delayed partial differential equation 



dtui 

dtU2 

dt(pR 

dt(t>o 
dtz 



CV • (yiVcjjR.) + ki 



0{yi) ui 



aiuf" — /3iMi 

a2U?''' - I32U2 

aR y2,tR + KflA((/)^«) - fca 3^1^ e{yi) ui - k^<pR<po 
ao y2,to + i^oM(l>'o) - k3(l>R<i>o 
-/i yi + h 2/2- 



(3.7) 



Recall that yi 



are the active cells and y2,t^ = y2{x,t — t^). The mechanisms behind 



BMU initiation are still not fully understood and we do not attempt to model them explicitly. 
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Instead, we initiate the BMU manually by perturbing the following fixed point of (3.7 1 



ui{x,t) 

U2{x,t) 

<j)o{x,t) 
z{x, t) 



U2,ss 




100. 



(3.8) 



To initiate the BMU we proceed now as follows. We leave the osteoclast field at steady-state 
ui^ss everywhere except for a confined region U where we add a few active cells ui^pertix) > for 
X & U. We assume that there are initially no active osteoblasts and that their density equals U2.ss 
everywhere. This is consistent with the assumption of the bone marrow being a precursor reservoir. 
The initial RANKL field is of great importance for the model because it is responsible for both 
targeted and strain-derived steering of the BMU. In fact, since neither the strain fields nor the 
osteocytes (which are responsible for the mechanotransduction by means of RANKL expression) 
are modeled explicitly as state-variables, possible damage sites and the principal stress directions 
have to be included in form of local perturbations of (f>R^pertix). Finally, we assume that there is no 
OPG present in the initial system and that the bone density is at 100% . In summary, the initial 
conditions are given by: 



Ul{x, t : 


= 0) 




U2{X, t : 


= 0) 




<f)R{x,t 


-0) 


= 4>R,pert{x) 


(j)o{x,t 


= 0) 


= 


z{x, t ~ 


0) 


100 



x ^n. 



(3.9) 



Since bone remodeling is a local process, we choose the domain large enough to avoid interactions 
of the BMU with the boundary. Note that for the BMU life spans considered hereinafter, large 
enough means at least one order of magnitude longer than the cutting cone. The corresponding 



Dirichlet boundary conditions for (3.71 are given in (3.8) with x € dVL 



Three comments regarding equations (3.7) - (3.9) are in order. First, we draw attention to the 



fact that the osteoblast and the bone density equations are ordinary differential equations and can 
be integrated explicitly. In particular, for the U2 equation we get 



U2{x,t)^U2.,sse~^'' + a2 / eP^'^'-'^uf'{x,s)ds. 



(3.10) 



Second, the Heaviside function introduces a discontinuity into the equations, raising questions about 
the well-posedness of the PDF. It can be seen that the point (ui^ss, U2,ss, 0, 0) is not a stable fixed 
point of the system. In the situations of interest, however, j/i cannot be zero unless 4>r is as well, 
since the active osteoclasts are only present in the cutting cone. Hence, we do not encounter issues 
of non- uniqueness. The questions of uniqueness and stability of the PDF system for the general 
situation are of interest, and are the subject of current work. 

Third, we expect the osteoclast field ui and the RANKL field (pn to inherit the singular behavior 



of the Heaviside function in (3.7). In addition, the RANKL field also suffers from porous diffusion 
effects, which themselves are known to exhibit singular behavior. If the initial RANKL field is 
compactly supported in a region with a smooth boundary, this free surface may develop local 
corners and cusps in the course of the simulation llj. Indeed, if we allow to become negative 
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(dropping below some threshold), very little can be said about the regularity of the ensuing PDE. 
This is an interesting question in its own right and will affect how computations may be performed. 
However, at this present juncture, we restrict ourselves to non-negative RANKL fields. 



4. The ID model. Due to the complexity of the model and the multitude of unknown pa- 
rameters, we look at the ID version of (3.7) - (3.9 1 before proceeding to the computationally more 
expensive 2D case. Note that in one dimension (n = 1), the differential operators simplify as 
V I— > 9j; and A i— ^ dxx- Before solving the system numerically, we first use some ideas of scaling 
theory to get a better understanding of physiological and pathological remodeling regimes as well 
as the corresponding parameter sets. 



4.1. Parameter estimation and sensitivity analysis. The primary goal after having es- 
tablished the model (3.7) - (3.9 1 is to identify a - not necessarily unique - set of parameters that 
corresponds to a physiological remodeling regime. Once this is achieved, various combinations of 
parameters can then be modified to study the emergence of pathologies. Ideally, the physiological 
parameter set could be estimated on the basis of experimental data. However, since almost none of 
our 23 parameters can be matched with experimental findings, we are forced to adopt a different 
strategy. First, we consider the purely temporal model (3.2 1 and follow the reasoning in [53] to 
obtain meaningful values. In particular, the values for /3i can be estimated from experimental find- 
ings about the corresponding life spans of bone cells. Also, it is shown that the value of gi2 leads 
to unstable results outside of the interval [0.1,4] and that gn determines the overall dynamics of 
the cell populations. These facts, together with an estimation of the time delays {tn, to) PI IMIIT^ 
and the aim of having a ratio of M2,ss/'"i,ss ~ 100 [IS], lead to the choice of a^, /3i, gij, tn and 
to found in (B.l). The remaining parameters cannot be matched with experimental data and we 
determine their physiological values a posteriori. More precisely, we fix the parameters in (B.l I, 
run simulations (as described in Section 4.2 ) and vary the remaining unknown parameters until 
the following two criteria are matched: first, the numerical solution has to coincide spatially and 
temporally with the global dynamics of in vivo observations and second, the cutting cone has to 
stay compact and move at a fairly constant speed. The outcome of this approach leads to the values 
summarized in (B.2). 

Now that the physiological set is determined, we can investigate the sensitivity of the model to 
parameter changes. To alleviate this task, we decide to focus on pathologies in the RANK/RANKL/ 
OPG pathway only. In other words, we consider the (B.l) parameters from now on as fixed pa- 
rameters and merely consider variations in the remaining free parameters of (B.2). However, a 
systematic sensitivity analysis of the 13 free parameters is still a rather unrealistic undertaking. 
Instead, we employ a scaling approach to analyze which parameters are able to destabilize the 
physiological regime and lead to the emergence of pathologies. The essence of scaling theory is to 
non-dimensionalize the equations by finding well-chosen scales for all the state-variables as well as 
the time and space variables. This leads to scaled equations where each term decomposes into the 
product of a dimensional coefficient representing the term's magnitude and a dimensionless factor 
of order of unity. Once this is achieved, it is possible to rewrite the equation in a dimensionless 
form where all the non-dimensional factors are now preceded by so-called dimensionless groups that 
contain all the information about the terms' magnitudes. The dimensionless ID version of (3.7) 
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reads 

diiii = Giuf" - G2U1 - G3a{yidis4>R) - G3bids:Uidi4>R) + Gij^^ ^{vi)ui 

df^R = Gry,^i^+Gsd,^{^'jf)~G,jij^e{m)u,-GioMo (4-1) 

di4>o = Giiy2:to+Gi2diS:{(j)'}^) - Gi3^R(j)o 
diz = -Guyi + Gi5 y2- 

The dimensionless groups Gi and the corresponding scales can be found in Appendix [A] Note that 
all the state variables Mi, 0^, z as well as x and t are now dimensionless and we can directly compare 
the various terms to determine their relative importance. In other words, we are now able to look 
for the dimensionless groups and parameters whose perturbations have a big impact on the model's 
regime. 

From a biological point of view, the most significant quantity is the bone mass density z{x,t). 
It contains the key information about the outcome of the remodeling process, i.e. it determines 
whether we have excessive, normal or insufficient remodeling of the bone tissue. Since the outcome 
of the bone mass balance is determined by the activities of osteoclasts and osteoblasts respectively, 
we have to focus primarily on the dynamics of iii and U2- However, bearing in mind that U2 only 
depends on ui and that the fixed parameters are kept at physiological values, we are assured that 
the osteoclasts will recruit enough osteoblasts to replace the resorbed tissue. In other words, the 
key players in the remodeling process are the osteoclasts and at this point, we do not have to worry 
about the osteoblasts. The only restriction to bear in mind is that the number of cells admissible 
per area is limited due to the cells' finite sizes; we ensure this by only considering parameter ranges 
that respect the spatial limitation. Osteoclasts are governed by the competition of G3 (magnitude 
of migration) and G4 x+^p^ (magnitude of stimulation by RANKL) and we define their ratio as 
(refer to Appendix [A] for the scales) 

G3 , A \ Ci"iA$fi , A \ CYi 

i 1 :— 




G4 V ^RJ hUiLl V '^rJ hUiLl 



Physiological remodeling only occurs if the two terms are well-balanced, Fi sa 1. A first pathological 
scenario corresponds to Fi ^ 1, i.e. the BMU moves much faster than it can nourish its population 
and dies out. On the other hand, if Fi <^ 1, we have too many osteoclasts produced in the cutting 
cone and hence too many osteoblasts recruited in the back of the BMU. Depending on the RANKL 
and OPG production rates, this can lead to an excessive production of RANKL which in turn 
creates more osteoclasts etc. This positive feedback loop in the closing zone can be investigated 
by means of the (j)R equation. A poor balance of RANKL production and its inhibition by OPG 
can lead to the described dysfunction in the closing cone zone of the BMU. More precisely, we are 
interested in the ratio of the production of RANKL by osteoblasts (G7) and its inhibition by OPG 
binding (Gio): 

Gj aRY2 aRl32 . , 

Lrio k^^R^o aoks^R 

A high ratio F2 ^ 1 leads eventually to a singular behavior of the model (blow-up of the cell 
populations). Yet another pathological mechanism involves the OPG field in the closing zone and 
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can lead to an early termination of the BMU. More precisely, if we have high production of OPG 
(Gil) in combination with low RANKL inhibition (G13), i.e. if 

p _ Gil _ aoY2 _ P2 
Gi3 fcg^o^fl fca^fl 

is very big, 3> 1, then the OPG field can possibly outrun the cutting cone and inhibit the 
RANKL field ahead of the BMU. The resulting lack of stimulation for the osteoclasts of the cutting 
cone can then lead to the extinction of the BMU. Obviously, this phenomenon only occurs if the 
diffusion is high relatively to the BMU speed. 

4.2. Numerical Experiments in ID. Note that the cutting cone of resorbing osteoclasts 
stays well-confined during the whole remodeling process and the BMU remodels a length of ap- 
proximately 5 mm in 6.5 months. Therefore, the simulation satisfies our criteria for a physiological 
regime and validates the choice of parameters. Calculating the ratios defined in Section |4!T] we get 
Fi — 0.83, r2 = 1.1 • 10^"^ and Fa = 2.7 • 10^'^. This is consistent with the previous discussion of 
parameter sensitivity. Indeed, Fi « 1 corresponds to a well regulated resorption activity, F2 ^ 1 
indicates a well-balanced RANKL distribution in the closing zone which is necessary for a confined 
cutting cone, and F3 <C 1 confirms that there is no risk of early termination due to excessive OPG 
production and diffusion. Finally, we point out that the scale estimations in Appendix A are in 
agreement with the simulation in Figure ??. 

Using the same set of physiological parameters, we investigate now the situation where a BMU 
starts off in the middle of two zones of high RANKL concentration (this corresponds e.g. to the 



situation of two adjacent micro- fractures) . Figure 4.1 illustrates how the cutting cone splits into 



two parts and remodels each zone separately. In particular, the BMU remodeling the higher peak 




Fig. 4.1. Physiological remodeling II. OC=osteoclasts, OB=osteoblasts, Z=bone mass. The length of the 
domain is 15 mm and the OB scale is to be multiplied by 10*. Note that the remodeling mechanism is adaptive: the 
higher RANKL peak at t = leads to more remodeling, see Z at t = 200. Parameter set and corresponding Ti as in 
Figure ?? . 
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is more active as can be seen in the bone density evolution. In other words, the rcmodcUng is 
adaptive: the bigger the damage and hence the RANKL expression, the higher the turnover in 
bone tissue. 

The remainder of this section is dedicated to pathologies. A first type of BMU malfunctioning 
is excessive bone remodeling and can be induced by two different imbalances. If we decrease the 




Fig. 4.2. Excessive remodeling. 0C= osteoclasts; length of the domain is 10 mm. A Increased osteoclast 
recruitment and lower RANK-RANKL binding saturation lead to a larger but compact cutting cone i n a s table 
reqirn e. The ratios are V\ = 2.6- 10~^, r2 = 1.3- 10~^ and Fa = 3.2- 10~^. The parameter set is given in \B.l^ and 
except for A = 2 and fci = 9 ■ 10~^. B Very low OPG production by osteoblasts in the closing zone lead to a 
slow and unconfined cutting cone. Positive feedback leads to instabilit y in the clo sing zone. The ratios are Fi = 52.2, 
F2 = 1.3 ■ 10"'^ and F3 = 1.7 ■ 10""^. The parameter set is given in SBTw and j^B.gp except for ap = 2 ■ 10~® after 
i = 60 days (ao is kept high in the beginning to avoid numerical instabilities in the initiation zone). 



ratio of osteoclast migration versus stimulation, i.e. if we choose the free parameters such that 
Fi <C 1, then more osteoclasts and hence osteoblasts are recruited and therefore the amount of old 
bone tissue that gets remodeled is expected to be much higher. If we simultaneously ensure that the 
feedback loop parameter is small, r2 <i; 1, we can avoid instabilities in the closing zone and expect 
an overall stable regime. These predictions are confirmed in the experiment illustrated in Figure 
|4.2| A. Note in particular that the cutting cone, even though much longer, stays confined and no 
instabilities occur. However, instabilities can no longer be avoided if excessive remodeling is caused 
by unbalanced RANKL/OPG production in the closing zone. In order to illustrate this, we pick a 
parameter set such that Fi « 1 but F2 ^ 1. As shown in Figure 4.2 B, the cutting cone is normal, 
but the excessive RANKL production in the closing zone leads to recruitment of a new generation 
of osteoclasts behind the cutting cone. These osteoclasts attract in turn more osteoblasts which 
produce more RANKL, and the resulting positive feedback loop leads to well- visible instabilities. 

Yet another pathological scenario is the early termination of the remodeling process, i.e. the 
extinction of the BMU before its mission is accomplished. Here too, we distinguish two different 
causes. If we choose Fi ^ 1, then according to our discussion in Section [41] the osteoclast popu- 
lation will die out due to deficient stimulation. Consequently, the whole BMU slowly disappears. 
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t=Odays t=15days t=35 days 




Fig. 4.3. Insufficient remodeling. OC=osteoclasts; length of the domain is 5 mm, A Decreased osteoclast 
recruitment and higher RANK-RANKL binding saturation lead to a vanis hing cuttin g con e. The ratios are Fi = 
10.5, F2 = 7.2 ■ 10"'* and Fa = 1.8 • 10""^. The parameter set is given in _?[ ) and except for A = 20 and 

fcl = 3 • 10^^. B High production and diffusion of OPG leads to annihilation of the RANKL ahead of the BMU 
and lack of stimulation le ads t o BM U ex tinction. The ratios are Fi = 0.58, F2 = 1.1 ■ 10"* and F3 = 10.5. The 
parameter set is given in \B. l\ and \B.l^ except for ap = 1 ■ lO"'^, ks = 1.5 • 10~^ and an = 10~^. 



see Figure [473| A. But early termination is also possible if osteoclasts respond well to RANKL stim- 
ulation: if the OPG production by osteoblasts largely exceeds the RANKL expression (Fa ^1) 
and if the OPG diffusion is very high, then the excess of fast spreading OPG reaches the RANKL 
ahead of the cutting cone and annihilates the osteoclast stimulation. Figure |4.3| B illustrates how 
the resulting lack in BMU stimulation can lead to early termination of the remodeling process. 

5. The 2D Model. We extend the model now to two space dimensions to gain a better insight 
into the dynamics of trabecular remodeling. Let € denote a rectangular domain representing 
the surface of a flat trabecula. The local cell densities of osteoclasts ui and osteoblasts U2 are 
denoted by Ui{x,t), where x = {x,y) G fi, i = 1,2. The RANKL field is denoted by 4>u{x,t) and 



the OPG field by <j)o{x,t). The governing equations are still given by (3.7) - (3.9) and V and A 
are now the Divergence and Laplace operators in 2D. Since the width of a trabecula is small in 
comparison to its length ii28| , and since the bone tissue is separated from the bone marrow through a 
canopy of bone lining cells ^21, vertical losses of RANKL and OPG are negligible (see also section|3|. 
This then justifies the use of a two dimensional diffusion equation to model the spread of chemicals 
across the surface of the trabecula. Note that we use the nonlinear, porous version of diffusion 
because the trabecular surface is very irregular and diffusing chemicals constantly interact with the 
components of the bone matrix as well as adjacent bone cells. In the remainder of this section, 
we present two numerical experiments on trabecular remodeling in a physiological regime. The 
first experiment is a regular micro-fracture remodeling and the second one illustrates that OPG 
plays an important role of counterbalancing the effects of RANKL and hence as a regulator for 
BMU-internal cell dynamics. More 2D-experiments in both physiological and pathological regimes 
together with a more detailed biological analysis of the results can be found in the accompanying 
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article [35]. Note finally that even though the scaling approach adopted for the ID case in Section 



4.1 looses its general validity in 2D, it can still be used to narrow down the plausible parameter 
ranges. 

5.1. Numerical Experiments in 2D. The following experiments are based on the model 



( 3.7 ) - ( 3.9 1, only the time delay terms in the RANKL and OPG equations are replaced by a^j y2,t^ ^ 



y2{x^E{x,t,t^), where 

, , \ _ ( ^ if y2ix, S) > for some 5 e [0, t — t^] andt > , , 

^[x,t,t^) - I Q otherwise. ^^'^^ 

This means that if at a certain location there exists an active osteoblast older than t^;, then all 
the active osteoblasts at the same location produce the respective chemical, independent of their 
age. This particular source term is practically useful because it does not require the use of delay 
differential equation solvers and hence improves both the computational cost and the stability of 
the algorithm. Furthermore, it is a reasonable approximation to the original version of the delay 



term y2,t^ as shown in Figure 5.1 In fact, considering the passage time of the cutting cone in the 
case of a physiological ID experiment shows that the latter is very short relatively to the time scale 
of the osteoblast dynamics. In other words, it is reasonable to assume that all the active osteoblasts 
at a specific location are of roughly the same age. In addition, the delay times are such that 
g-/32tij « 1, and we conclude that y2{x,t) 5(a;,t, t^) is indeed a reasonable approximation for ?/2,t„- 
All the simulations are performed in Matlab by means of a second order finite difference scheme in 




Fig. 5.1. Normalized population dynamics of active osteoclasts (dashed line) and active osteoblasts (solid line) 
at X = 3.7mm in the experiment Physiological Remodeling I (see Fig. ??). We see that the passage time of the 
cutting cone is very short relatively to the time scale of the osteoblast dynamics. Therefore, one can assume that all 
the active osteoblasts are of approximatively the same age. 



space and the built-in solver ode45 in time. 
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First, we demonstrate the effect of RANKL on BMU steering along a micro-fracture. The me- 
chanically damaged osteocytes adjacent to the fracture create a path of membrane-bound RANKL 
as depicted in Figure [5^ at time t = 0. In the course of the simulation, the RANKL-guided BMU 





Fig. 5.2. RANKL field simulating microfracture in trabecular bone. Prior to the simulation, damaged 
osteocytes along the fracture express membrane-bound RANKL leading to the initial field at t=0 days (white indicates 
high concentration). At time t=160 days, the RANKL field has almost entirely disappeared after having bound to 
both OPG and RANK receptors on osteoclasts. Since RANKL is membrane-bound, the diffusion is very low. Length 
of domain is 3 mm. 



remodels the fracture and the RANKL disappears due to RANK-RANKL binding, leading to the 
final snapshot after t = 160 days. We initiate the BMU by introducing a confined aggregation of 
active osteoclasts at the bottom of the microfracture at time t = 0. The first panel in Figure 5.3 
shows the subsequent motion of the cutting cone: the bright area represents the region of active 
osteoclasts which move towards the top of the fracture. The osteoblast dynamics are depicted in 
the second panel: osteoblasts are recruited by active osteoclasts and produce new bone matrix in 
the areas where the cutting cone has already resorbed the bone. The third panel shows the OPG 
field: it is produced by mature osteoblasts and hence lags the cutting cone. In particular, OPG 
is not membrane-bound and diffuses across the trabecula at a fairly high speed. The last panel 
depicts the evolution of the bone density. Note finally that the cutting cone stays well-confined and 
the BMU moves at constant speed over 2 mm in 5 months. This is in agreement with experimental 
observations [55] and thus provides a validation of the chosen physiological parameter set. 

The second experiment is an extension of the ID experiment on the possibility of BMU branch- 
ing in the case of multiple micro-fractures (see Figure 4.1 1 . More precisely, we want to find out if 



a BMU can split into two separate BMUs and if it can, then we want to investigate the existence 
of preferential branching directions. We start off with the initial RANKL field from the previous 
experiment and add a secondary branch which deviates by 45° from the primary branch as shown 
in the snapshot RANKL fwd at t = 15 in Figure 5.4 Again, an initial aggregation of osteoclasts is 
placed at the bottom of the microfracture and as time progresses, this cutting cone moves towards 
the top of the fracture. Similarly to the ID experiment, the BMU splits into two individual parts 
which remodel both branches separately. Interestingly, if one repeats this experiment with the 
secondary branch deviating at 135° rather than 45° (see RANKL bwd), the BMU remodels the 
primary but not the secondary branch. In fact, the RANKL in the secondary branch is annihilated 
by OPG-RANKL binding before the cutting cone reaches the branching location, and in this way - 
due to the resulting lack of RANKL stimulation - the osteoclasts do not deviate from the primary 
branch (see OC bwd). In summary, the OPG production in the back of the BMU prevents the BMU 
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Fig. 5.3. Steering of BMU along microfracture. OC; aggregation of osteoclasts (cutting cone) moving from 
the bottom of the domain to the top along the RANKL gradient. OB; osteoblasts, rebuilding the bone in the wake 
of the cutting cone. OPG; diffusing OPG field. Z; evolution of the bone mass density. Outline of initial RANKL 
field (micro-fracture) is highlighted for reference; length of domain is 3 mm; black corresponds to low, white to high 
concentrations. 

from turning around and remodeling the previously remodeled tissue. This is in good agreement 
with experimental results obtained by means of microcomputed tomography (MCT) imaging by 
[lOj . For a more detailed discussion of this branching phenomenon we refer to |32j . 

6. Conclusion and Outlook. We have established a novel mathematical model of bone 
remodeling at cellular level. Based on a previous temporal model by Komarova et al. we developed 
step by step a spatio-temporal model describing both the osteoblast and osteoclast populations 
as well as the dynamics of the RANKL and OPG fields. The complete model has been shown to 
successfully recapitulate the overall dynamics of a single BMU as well as the distinct features of 
the cutting cone. Scaling was used to investigate the importance of the various model parameters 
and to motivate experiments on pathological remodeling. 

A strong feature of our model is the possibility to investigate the role of the spatial RANKL and 
OPG distributions in the osteoblast-derived paracrine control of osteoclasts. Even though there is 
a strong consensus in the experimental literature about the importance of the RANKL/OPG ratio, 
the following apparent inconsistency has to our knowledge not yet been addressed. In fact, the 
cutting and closing cones are spatially disconnected and hence osteoblasts appear when osteoclasts 
are already gone. So how can osteoblasts possibly play a key role in the osteoclast control ? Our 
results show that the spatially distinct distributions of the RANKL and OPG fields provide the 
missing link: by expressing the diffusing RANKL- inhibitor OPG, osteoblasts have an indirect means 
of control over the activity of osteoclasts and hence the extent of remodeling and the direction of 
movement of the whole BMU. The 2D version of the model is particularly suited to describe 
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RANKLtwd OCtwd RANKLbwd OPGbwd OCbwd 




Fig. 5.4. Forward versus backward branching. In addition to the primary micro-fracture, secondary fractures 
branching at 45° and 135° are added (see RANKL fwd and RANKL bwd at t = 15j. OC fwd.- Remodeling in the 
forward direction is successful, the BMU splits into two parts and remodels each branch separately. OC bwd; 
Remodeling in the backward direction is unsuccessful, the BMU only remodels the primary branch. The reason for 
this is annihilation of RANKL by OPG along the secondary branch before the cutting cone reaches the branching area 
(see RANKL bwd and OPG bwdj. The resulting lack of osteoclast stimulation prevents the BMU from branching 
and the backward branch are not remodeled. Length of domain is 3 mm. 



trabecular remodeling in the case where the local curvature of the trabeculae is negligible. Regarding 
cortical remodeling, it is likely that the 2D model provides a good qualitative approximation in the 
case when the BMU moves within the same plane of the cortical tissue. Nevertheless, in order to 
draw quantitative conclusions, a full 3D formulation together with a few modifications of the model 
assumptions become necessary. 

For future investigations, the model presented in this article provides a promising starting 
point. Besides an improvement of the numerical scheme and an extension to three dimensions for 
cortical remodeling, we also plan to improve our results by adding the natural decay rates as well 
as appropriate stochastic terms to the RANKL and OPG equations. In fact, since the production 
of messenger molecules by cells are subject to fluctuations, the use of noisy RANKL and OPG 
sources is expected to improve the model predictions in view of the often very irregular BMU 
evolutions observed in vivo. Further model improvements might be achieved by describing precursor 
cells as independent state variables and by including other important regulating factors such as 
Sclerostin, TGF-/3 and PTH as state variables. However, the resulting increase in complexity would 
further compromise the balance between reliability and realism: the parameter-fitting for the current 
model already presents a substantial challenge and the addition of more unknown parameters would 
certainly not improve the model's quantitative reliability. Regarding the mechanical factors, model 
improvements can be achieved by taking into account the local curvature and by coupling the model 
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to existing finite element models describing the elastic properties of the tissue. 



Appendix A. Dimensionless groups and scale estimations. Due to the multiple time 



and length scales in the system (3.7 1 as well as the occurrence of two different zones (the cutting 
zone and the closing zone), we have to abandon the idea of finding a consistent non-dimensional 
version of the original equation with a single set of scales. However, we can still transform (3.7 1 into 



a dimensionless equation where the dimensionless factors preceded by the dimensionless groups are 
all of the order of unity - all we have to do is to use different scales for different terms. Furthermore, 



the structure of the resulting equations (4.1 ) with respect to the time derivative implies that even if 



we cannot identify a single time scale for each equation separately, we can still compare the terms 
on the right hand side because the ratios of the form Gi/Gj are independent of the time scale. The 



dimensionless groups in (4.1) are defined as 



Gi=T (aiC/f"~^ 
G4 = T fci 



Gi2 = T 



G2 = r/32 

Q — 2^ k2Ul 

Gi^=Th^R 



G6^TI32 

Gio -T/c3$o 



G?.b-1 \u^LnL^) 

G-! = T 
Gil 



$0 



(A.l) 



Gi5 = r^ 



Except for i £ {8, 12, 15, 16}, most of the terms Gi play a significant role in only one of the two 
remodeling domains: i £ {1, 2, 3, 4, 5, 6, 9} in the cutting cone and i G {7, 10, 11, 13} in the closing 
cone. This has to be kept in mind when looking for the correct scales. The capital letters Ui, Yi, 
<i>fl, $0 and Z are the scales of the corresponding state variables. Li, Ljj, Lq are the length scales 
of the osteoclast, the RANKL and the OPG fields respectively, and 'tn scale the RANKL field 
at the tip of the cutting and in the back of the closing zone respectively. A^/j is the difference of 
the RANKL concentration between the front and the back of the cutting zone. 

In the remainder of this section, the various scales are now briefiy justified. Since we assume 
physiological remodeling conditions, the length of the cutting cone and the number of its constituent 
active cells are supposed to be preserved. Therefore, the length scale Li equals the initial length of 
the cutting cone and Yi « ui^pert with ui^p^rt the initial perturbation added to the steady-state pool 
of passive cells at time t — 0. Consequently, Ui ~ ui_ss + ui^pert- For the RANKL field we note first 
that the biggest change in concentration occurs in the cutting zone and hence the corresponding 
length scale is L/j « Li. Since physiological remodeling excludes excessive RANKL production by 
osteoblasts, the scale is dictated by the initial conditions, $ij w maXx£n\4'R{Xjt = 0)|. Next we 

estimate the passage time Tp it takes the cutting cone to move across its own span: Li is divided 

^2 

by the velocity of the osteoclasts to get Tp — ^^^^ . This expression allows us to eliminate Tp in 
the estimation for A$j^ 
field, 



Tpk2Ui xff^ ^-iid we obtain, respecting the positivity requirement of the 



A$J^ w min<^ 



C(l 



(A.2) 



The remaining RANKL scale is given by ~ — A$/{ + aji{to — t/s)i^2- Using the time 



Tp we get then for active osteoblasts Y2 ~ Tpa2Yi = 



Yi and hence U2 



"2„ 
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Since the OPG field is generated by active osteoblasts with a life span of 1//32, we get the estimates 
$0 ~ 0,0^2/ h and Lo ~ ^^^^f ■ Finally, the bone mass is scaled with respect to Z ^ 100. 

Appendix B. Parameter Values for ID Experiments. As discussed in Section 4.1 



we distinguish between fixed and free parameters. The former are unchanged throughout all the 
experiments and their numerical values are summarized in the following table. 



ai = 9.49 day-imm-°^ Q!2 = 4 day-^ f3i = 0.2 day-^ = 0.02 day" 

511 = 0.5 512 = 1 tii^A day to day 



(B.l) 



The free parameters are changed from experiment to experiment. For the physiological experiments 
we use the following set. 



^ = 7 . 10 mm^ mol" ^ day" 
an — 6 ■ 10~^ molday"^ 

eji = 2.5 eo = 1 



fcl = 3 • 10 ^ day-^ A = 5 molmm"^ 

ao ~ 1.5 • 10"'' molday"^ Kfl = 3.16 • 10^^ mm'ri + ^ mol^^'R day"^ 
= 1 • 10 molday"^ — 1.2 mm mol" ^ day " ^ 



KO = 10" 



•O+i mol^-'O day-i /i — 0.3 gday" 



/2 = 1.6-10-^ gday- 



(B.2) 



Appendix C. Parameter Values for 2D Experiments. 



511 



30 day- 

= 0.5 



— 10 mm^ mol ^ day 
ao = 3 • 10^^ mol day" ^ 



510- 



mm mol day 



Q;2 = 4 day"^ 
512 = 1 

fci = 2.8 • 10-3 day-i 
Kji ~ 10^^ mm^Ti mol^ 
Kq = 10^^ mm^'O mol^ 



day" 
'O day- 



/?! = 0.1 day-i 

tn. — 5 day 

A — 50 molmm"^ 

fi? = 3 eo = 1 

/l =0.24gday-i 



/?2 = 0.02 day-1 

to = 15 day 

flfl = 10^^ molday"^ 

k2 = 4.6 • 10^^ mol day" 



/2 = 1.7 • 10- 



5 day 



(C.l) 
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